Molecular dynamics (MD) simulations provide invaluable insight into the properties of biomolecular and other complex systems such as macromolecular systems. Some problems that are beyond experiment can only be tackled by simulation. Unfortunately, accurate simulations of macromolecules are often computationally demanding and in many cases are not feasible due to the large number of particles that involve complex and long-range interactions. New approaches and improved computing platforms are needed to enhance the reliability of macromolecular simulations and to approach realistic time frames.
Biological and some other complex systems are multi-scale in nature. For example, the dynamics of proteins contain motions over different time scales, from atomic vibrations in the order of femtoseconds to collective motions at millisecond scales. FIG. 1 depicts the dynamics of molecules such as protein molecules, to illustrate the variation in time scales.
Traditional time stepping integrators (e.g. Verlet) commonly used in molecular dynamics (MD) are not able to address this time scale problem. A typical time-step for these methods is 1 femtosecond. This makes atomistic simulation of biomolecules computationally extremely expensive. Multi-scale numerical methods, in which the presence of fast scales does not affect the time integration of slow scales, are urgently needed for efficient simulation of large biomolecular systems. Such approaches, in theory, can essentially enhance performance of molecular simulation since the most computationally expensive long-range electrostatic interactions contribute to the dynamics on relatively long time scales (compared with internal vibrations) and thus ideally do not need to be calculated frequently. Also, such approaches enhance the data locality which makes them better suited for implementation on parallel computers than traditional MD schemes.
Multiple-time-stepping (MTS) methods are among the most popular methods of this type. In MTS methods, savings of computational time can be realized if the slowly varying forces due to distant interactions are held constant over longer intervals than the more rapidly varying short-range forces. Standard integration procedures in MD can then be modified by evaluating the long-range forces less often than the short-range terms. The ratio between frequencies of evaluation of the long-range forces (outer step-size) and short-range forces (inner step-size) measures the gain in simulation time and will be further referred as “the step-size ratio”.
FIG. 2 is a diagramatic explanation of the MTS idea, which splits the forces in a system into bonded “fast” forces and long range non bonded “slow” forces (which tend to be non-linear), evaluating the slow forces less frequently. For this, multiple timestepping integrators are required to solve modified ODEs (ordinary differential equations).
For biomolecular applications, the computational complexity of the fast, short-range force evaluations scales linearly in the number of atoms in the system, N, while it scales quadratically in N for the slow, long-range force evaluations. Furthermore, while the short-range fast forces are easy to compute in parallel, long-range slow forces require global data communication and hence are more difficult to parallelize efficiently. Therefore, in theory, MTS methods can dramatically speed up MD simulations by reducing the number of expensive slow force evaluations.
In practice, however, MTS methods such as the popular Verlet-I/r-RESPA (Grubmüller et al., 1991; Humphreys et al., 1994) suffer from severe resonance instabilities that limit practical performance gain (Ma et al., 2003; Izaguirre et al., 2001). For solvated biomolecular systems, one, for example, has to update expensive slow forces only four times less often than calculating cheap fast forces.
Another problem is that common time-stepping methods do not exactly sample from the target temperature even if the simulations are stable and are subject to a thermostat (Pastor et al., 1988; Bond and Leimkuhler, 2007). The resultant error can be controlled by increasing the frequency of updating slow forces and, in fact, calculating them much more often than required by stability considerations alone. This obviously reduces computational efficiency.
These resonance induced instabilities of impulse MTS methods have been reduced, through the introduction of mollified MTS methods by Izaguirre et al., 1999, giving improved linear stability by defining a slow part of potential energy at a time-averaged position. Further improvements have been achieved by weak coupling to a stochastic heat bath (Langevin dynamics) (Izaguirre et al., 2001). However, accurate simulations still put limits on the step-size ratio which to weaken non-linear instabilities (slow forces) must be in the range of 6-12 for solvated biomolecular systems (Izaguirre et al., 1999, 2001). Moreover, these improved methods can suffer from not reproducing system properties accurately and it may still be necessary to chose simulation parameters very carefully to provide a stable simulation.
It is desirable to provide a method and apparatus for simulation which overcome or at least mitigate some of the disadvantages of the prior art.